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SUMMARY 


The  objective  of  this  program  is  to  investigate  surface-patch 
modeling  and  to  develop  a  user-oriented  computer  algorithm  for  analyzing 
general  three-dimensional  radiation  and  scattering  problems  with  arbitrary 
surface  geometries  for  airborne /ground-based  applications.  The  existing 
surface-patch  algorithm  was  improved  and  tested  on  prolate  spheroids, 
circular  and  rectangular  cylinders,  and  boxes.  Good  numerical  results 
were  obtained  even  for  geometries  with  sharp  edges.  The  merits  of  the 
surface-patch  approach  in  solving  electromagnetic  problems  of  arbitrary 
geometries  have  been  demonstrated  for  conducting  bodies  whose  dimensions 
are  not  electrically  large.  Its  advantages  as  well  as  disadvantages  in 
comparison  with  the  wire-grid  approach  were  also  discussed.  Development 
of  a  new  surface-patch  algorithm  was  started  with  a  reaction  integral 
equation.  This  approach  required  intolerably  large  computer  time  in 
the  Integration  process  involved  in  the  matrix  computation.  The 
electric  field  integral  equation  was  then  used  to  overcome  this  problem. 
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ACCURACY  OF  THE  64-POINT  NUMERICAL  INTEGRATION 
ALGORITHM . 


SECTION  1 
INTRODUCTION 


There  are  presently  two  types  of  approaches  that  have  been  employed 
to  treat  the  scattering  and  radiation  of  arbitrarily-shaped  conducting 
bodies  that  are  not  electrically  large.  The  first  is  to  simulate  the 
conducting  surface  with  thin  wire  grids  and  the  second  is  to  approximate  the 
surface  with  patches.  Research  in  the  wire-grid  approach  has  been  very 
active,  yet  comparatively  little  effort  has  been  devoted  to  the  surface 
patch  approach.  Surface-patch  modeling  of  arbitrarily-shaped  dies  has  been 
investigated  by  Oshiro,  et  al.  (1,2],  Albertsen,  et  al  [3’  j  Wang 
[4]. 

The  objectives  of  this  research  program  are  to  investigate  surface- 
patch  modeling  and  to  develop  user-oriented  computer  programs  for 
analyzing  general  three-dimensional  radiation  and  scattering  problems 
with  arbitrary  surface  geometries  for  airborne  and  ground-based 
applications.  The  tasks  include  the  evaluation  and  improvement  of  the 
existing  Georgia  Tech  algorithm  and  the  development  of  computer  programs 
for  thin-shell  and  antenna  problems.  This  interim  report  summarizes  the 
accomplishments  during  the  first  year  of  this  two-year  research  program. 

The  algorithm  improvement  was  primarily  directed  toward  the 
reduction  of  execution  time  and  the  amount  of  central  memory  required  in 
the  numerical  computation.  Numerical  testing  was  carried  out  for 
prolate  spheroids,  circular  cylinders  and  rectangular  cylinders  (boxes) 
in  order  to  determine  the  accuracy  of  the  algorithm.  Also,  a  computer 
program  employing  the  electric  field  integral  equation  was  developed  to 
treat  thin-shell  problems  which  can  not  be  computed  using  the  existing 
Georgia  Tech  algorithm. 
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SECTION  II 

EVALUATION  AND  IMPROVEMENT  OF  EXISTING  ALGORITHM 

A.  Basic  Formulation 

The  existing  computer  program  at  Georgia  Tech,  as  discussed  in 
Reference  4,  employs  a  magnetic  field  integral  equation  (MFIE)  as 
follows 


J 

— s 


(r) 


2~~  n  X  x  (r,r’)ds’ 

2n  x  H  (_r)  for  r;eS  ( 


(1) 


where  *  surface  current  density  on  the  surface  S, 
n  =  an  outward  unit  vector  normal  to  S, 

*  [exp(- jk)  jr  -  r'  ))]/  J  jr  -  _r'  j , 

jr,  _r*  »  vectors  from  the  origin  to  the  field  and  source  points, 
inc 

H_  *  incident  magnetic  field  intensity, 

J-  =  principal  value  integral,  and 
V'  =  gradient  operator  with  respect  to  the  primed  coordinates. 

For  the  convenience  of  the  present  discussion,  it  is  desirable  to 
review  the  matrix  generation  process  involved.  We  define  the  operator  L 
as 


L(J.s(r))  =  J.g(r) 


1 

2n 


fS 


J^r'lx  V'  cpds 1 


(2) 


We  now  let 


and 


J  (r)  =  I' LP  \r)  +  I  2P  2(r), 
n  n  n  —  n  n 


Pi(r)  -  U1 


for  reAS 
—  n 


*  0  elsewhere, 


1  2 

where  U  and  U  are  orthonormal  unit  vectors  on  AS  defined  by 
— n  — n  n  J 


(M 

(5) 


n  x  HinC(£) 

I  n  X  HinC(r) 


(6) 


U  2  =  n  x  U  1. 
— n  n  — n 


(7) 


Using  these  relations,  Equation  (1)  can  be  written  in  the  form 


L  [P^(r)l  -  P*(r)  -  ds', 

n  —  n  Lv  n  J 


AS 


and  we  define  the  weighting  function  as 


(8) 


W  i(r)  -  «(r  -  r  )P  i(r), 
— m  —  —  m  m  — 


(9) 


and  the  scalar  product  as 

<P,Q>  =  f  P  •  £  ds.  (10) 

*'s 

Thus,  by  applying  the  weighting  function  to  Equation  (8),  we  have 
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(11) 


(r),  L [P 1(£)]> 


n 

tn 


X  {U  jX  /  v'*(r  r')dsM  . 

-Ti  J  "in. 

Asn 

The  magnetic  field  integral  Equation  (1)  can  be  used  to  generate  a 
system  of  linear  equations  by  expanding  the  unknown  surface  current 
according  to  Equation  (3)  and  performing  the  scalar  product  according  to 
Equation  (11).  The  resulting  system  of  linear  equations  can  then  be 
denoted  by 

2  «  *  i  <  i 

y  v  I  jZ  1  j  =  V  1  (12) 

“  ■"  n  m  n  m 

j=l  n=l 

m  =  1} • • *N 
i  =1,2 

where 


Z  ±  3  =<W  i(r),L  [P  3(r)l> 
m  n  N— m  —  n 

Vi  =<W  i(£),  2n  X  HlnC(r)>  . 
m  -m 


To  solve  for  the  unknown  surface  current  I  one  must  first 

13  n 

compute  the  matrix  element  which  involves  an  integration  process 
according  to  Equation  (11).  The  existing  Georgia  Tech  computer 
algorithm  is  based  upon  the  formulation  which  is  summarized  in  the 
preceding  equations.  The  methods  of  evaluation  of  the  matrix  elements 
and  matrix  solution,  together  with  the  techniques  which  have  been 
developed  to  improve  the  computer  algoirthm  will  be  discussed  in  the 
following  paragraphs. 

B.  Matrix  Symmetry 

Execution  time  and  central  memory  size  are  primary  considerations 
in  the  present  analysis.  It  is  desirable  to  explore  the  possibility  of 
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attaining  a  symmetrical  property  for  the  square  matrix  ,  A 

mn 

symmetrical  matrix  requires  only  about  half  the  execution  time  and 
central  memory  as  does  an  unsymmetrical  matrix.  Thus,  symmetry  will 
make  it  possible  to  compute  the  EM  field  for  scattering  and  radiating 
bodies  of  larger  size. 

The  square  matrix  can  be  arranged  in  the  following  form 


[Z 


i  J 

m  n 


(Z 


1  1 
m  n 


] 


[Z 


2  1 
m  n 


] 


.12 

tZmn 


[Z 


2  2 
m  n 


(13) 


which  consists  of  four  sub-matrices  each  of  which  has  a  NxN  dimension. 
Letting 


<  -  m  +  (i-l)N, 
*  -  n  +  (j-l)N, 
we  have 


(14) 

(15) 


k  i  j 

C  ■  Z  J 
l  m  n 


(16) 


For  [Z^  ]  to  be  symmetrical,  it  is  necessary,  by  definition, 
mn 

that 


(17) 


or 


z  1  J  -  c*  -  c" 

m  n  £  x 

-  Z  ^  1  . 
n  m 


(18) 
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Substitution  of  Equation  (11)  into  Equation  (18)  yields 


n 

m 


x  (V* 


(r  ,ir '  )ds' } 


V 


X 


(r  ,r_')ds’} 


(19) 


which  is  the  condition  for  a  symmetrical  matrix  [Z^  ]  .  In  general, 

mn 

Equation  (19)  can  not  be  satisfied  and  therefore  symmetry  for  the  matrix 
can  not  be  readily  achieved. 

Galerkin's  method,  which  chooses  identical  sets  of  basis  and  weighting 
functions,  is  more  likely  to  yield  symmetrical  matrices.  In  the  present 
case,  this  can  be  explored  by  selecting  equal  patch  sizes  and  weighting 
functions  as  follows 


— m 


(r) 


(r). 


(20) 


The  condition  for  matrix  symmetry  in  this  case  is 


U 


i 


— m 


n  x 
in 


4^  *  /  f  V'<Kr,r-)ds'ds} 

> sm  Jl sn 

x  {^tn  X  I  {  v'4>(r,  r')dS'ds}, 

>  Sn  Xsm 


(21) 


which  can  be  reduced  to 


-m 


n  x  =  l4 
m  -n  -n 


n 

n 


U 

m 


(22) 
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or 


iP  •  U1+<51  (-I)"5 
— n  — m 


(23) 


Unfortunately,  Equation  (23)  can  not  be  satisfied  on  a  general  three- 
dimensional  surface.  If  we  let 

1^  =0  for  n  =  1,  ...  N  (24) 

n 

in  Equation  (4),  Equation  (23)  is  reduced  to 


This  is  similar  to  wire-grid  modeling,  in  which  the  direction  of  the 
current  is  along  the  pre-determined  wire  structure.  The  condition 
stated  by  Equation  (25)  is  still  too  restrictive  to  be  satisfied  on  a 
general  closed  surface. 

The  above  discussion  can  also  be  reviewed  on  the  basis  of  linear 
algebra,  which  attributes  the  symmetry  of  the  matrix  to  the  self¬ 
conjugate  property  of  the  linear  operator  L  [5],  As  long  as  one  selects 
an  orthonormal  basis  for  a  space,  then  self-conjugate  operators  are  in  a 
natural  one-to-one  correspondence  with  symmetrical  matrices.  The 
difficulty  in  attaining  matrix  symmetry  therefore  appears  to  be  due  to 
the  magnetic  field  integral  equation  for  which  the  operator  L  is  not 
self-conjugate.  It  is  therefore  necessary  to  modify  the  operator  L  ,  or 
the  magnetic  field  integral  equations,  before  a  further  step  to  attain 
matrix  symmetry  is  attempted. 

In  the  wire  scattering  and  radiation  analysis,  apparently  the  only 
algorithm  having  a  symmetrical  matrix  is  the  reaction  integral  equation 
developed  by  Richmond  [6],  The  original  operator  was  not  self-conjugate 
and  the  matrix  therefore  asymmetric.  However,  by  approximating  the 
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surface  integral  over  the  tubular  expansion  dipole  with  a  line  integral, 
the  matrix  element,  or  the  mutual  impedance  between  two  dipoles,  was 
made  symmetrical.  Other  computer  programs  do  not  have  a  symmetrical 
matrix  in  their  computation.  However,  matrix  reduction  for 
geometrically  symmetrical  scatterers  and  radiators  has  been  recognized 
and  practiced  such  as  in  theWAMP  program  17],  the  work  by  Sancer,  et  al. 
[8],  and  Tsai,  et  al.  (9],  In  the  present  algorithm,  matrix  reduction 
for  a  symmetrical  scatterer  was  implemented  and  will  be  discussed  next 
in  this  section. 


nmetrical  Scatterers 


When  the  direction  of  propagation  of  an  incident  plane  wave  is  in 
the  plane  of  symmetry  of  a  conducting  scatterer,  it  is  recognized  that 
some  symmetric  behavior  must  exist  in  the  induced  current  on  the  surface 
of  the  scatterer.  This  property  has  been  taken  advantage  of  by  various 
authors  in  their  effort  to  reduce  the  central  memory  size  required  in 
the  computer  run  [7-9]. 

Without  loss  of  generality,  Cartesian  coordinates  can  be  set  up  so 
that  the  plane  of  symmetry  coincides  with  the  XZ  plane  as  shown  in 
Figure  1.  The  polarization  of  the  incident  wnve  is  assumed  to  be  either 
parallel  or  perpendicular  to  the  XZ  plane.  Arbitrary  polarizations  can 
be  decomposed  into  two  components,  one  parallel  and  the  other 
perpendicular  to  the  XZ  plane.  The  overall  scattering  problem  can  then 
be  treated  by  superposing  the  fields  due  to  these  two  component  incident 
fields. 

At  two  symmetrical  points  I  and  I+N/2  In  Figure  1,  the  components 
of  the  induced  surface  currents  exhibit  the  following  relationships 
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when  11.  is  perpendicular  to  the  plane  of  symmetry  (XZ  plane).  In 

1  ix  11 

Equation  (26),  J  denotes,  for  example,  the  x  component  of  Ir  P  in 

I  - n 

Equation  (4). 

When  is  parallel  to  the  plane  of  symmetry,  the  induced  current 
on  the  scatterer  has  the  following  property 
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lx 
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lz 

I 
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-Jlx 

I+N/2 

Jly 

I+N/2 

lz 


-J 


-J 


I+N/2 

2x 

I+N/2 

2y 


I+N/2 

J2z 

I+N/2 


(27) 


Equations  (26)  and  (27)  show  that  the  number  of  unknowns  in 
Equation  (12)  can  be  reduced  from  a  total  of  2N  to  N.  This  can  be 
carried  out  by  the  following  reduction  process 


r  j  7*j 


N  2  N/22  .  .  .  N 

£  £  ‘J  z»  *  £  £  li +  £  Vi- 

n*l  j-1  n=i  j=l  n=N/2  +  1  j=l 

N/2  2 


(28) 


n/ 1  £  .  r  J  •  <  I  1 

+  Z’m(n+N/2)j  , 


where 


Z’i(n+N/:>  =  (^>  L[^IP*  ^  (^)]>>  <29> 


where  is  a  dyad  defined  as 

=IP 


10 


cJT„  -  (-D(j-1){(-i)Ip+j+1«+  (-i)Ip+j+2^+ 


-ip 


(30) 


Ip  =  1  if  H  //  x-z  plane 
=  2  if  H*  J_  x-z  plane 


Equation  (12)  can  now  be  written  as 


N/2  2 

E  E 

n=l  j=l 


rJ 


[Z 


mn 


+  Z 


.ij 

m(n+N/2) 


i  -  1,2 

m  “  1 .... N/2 


(31) 


The  total  number  of  equations  is  then  reduced  to  one  half  of  the 
original  number  and  the  matrix  size  to  one  fourth  of  the  original  size. 

D.  Matrix  Computation 

The  matrix  element  2^  ,  according  to  Equations  (11)  and  (12), 

mn 

involves  an  integral  of  the  following  form 


k 


V'<f>  (r  ,  r/)  ds\ 


(32) 


Computation  of  the  integral  1^  is  very  time  consuming  if  carried  out 
numerically.  However,  it  has  been  observed  that 


1  ^  V'*  (r  ,  r  )  AS  , 

— n  "“TO  ii  n 


(33) 


where  is  the  position  vector  to  the  center  of  the  nth  patch  and  ASn 
is  its  area.  The  approximation  in  Equation  (33)  Improves  as  the  shape 
of  the  nth  patch  approaches  an  equiangular  triangle. 
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Figure  2  shows,  without  loss  of  generality,  the  nth  triangular 
patch  coincidng  with  the  X'-Y'  plane,  one  vertex  on  Y'  axis,  and  the 
opposite  side  on  X*  axis.  The  integrand  can  be  simplified  by  the 
following  approximation 


4>  (r  ,  r') 


-jkR 


mn 


[l-jk  (R-R  ) ]/R» 
mn 


(34) 


—  i|cr 

which  is  a  two-term  Taylor  expansion  of  e  in  the  neighborhood  of 

R  ,  where  R  ■  I  r  -  r'  I  ,  R  *  I  r  -  rlasin  Figure  3.  The  validity 
of  Equation  (34)  depends  on  the  following  condition 


(35) 


Since  the  patch  size  is  limited  to,  say,  X/3  or  less  on  each  side  of  the 
triangle,  it  is  generally  true  that 


k  |  R-R  |  <  1  .  (36) 

1  mn 

It  can  be  shown  that  the  self  patch  integral  I  °  vanishes  when  the 

— n 

patch  Is  equilateral  and  small  enough  to  satisfy  Equation  (35). 
Therefore 


n  t 
1=0 


(37) 


In  proving  Equation  (37),  it  was  noted  that  Equation  (35)  was  well 

satisfied  in  the  immediate  vicinity  of  R  ,  where  |R  -  R  I  **«  |  R  |  ~  0. 

mn  nm 

When  n--£n,  two  separate  patches  are  involved,  the  integral 
becomes  considerably  more  complicated  but  can  be  approximated  in  closed 
form  as  follows 


X  m  n. 
— n 


(1  +  jkR 


mn 


-jkR 


mn 


(x'l 


y’V’ 
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Figure  3.  Geometry  of  the  mth  and  the  nth  patches. 
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where  x,y,z  denote  x  ,  y  ,  z  ,  respectively. 

m  m  m 

Although  Equation  (38)  Is  an  accurate  expression  for  the 
computation  of  ,  it  is  quite  cumbersome.  In  practice,  the  simple 
expression  of  Equation  (33)  is  quite  satisfactory.  In  fact,  all  the 
numerical  data  included  in  this  report  are  based  on  Equation  (33).  The 
CPU  time  required  in  generating  the  matrix  elements  is  small,  being  on 
the  order  of  100  seconds  for  a  192x193  matrix  on  the  CL)C  CYBER  74 
computer  at  Georgia  Tech. 

E,  Interior  Resonance 

Interior  resonance  in  electromagnetic  scattering  has  been  observed 
in  two  types  of  problems,  one  due  to  the  simulation  by  wire-grids 
[10,11],  the  other  being  related  to  the  integral  equation  used.  In  both 
cases,  the  solution  to  an  integral  equation  is  not  unique  at  the 
resonance  frequency  of  the  interior  cavity.  Thus,  a  solution  to  the 
interior  problem  emerges  while  the  exterior  problem  is  being  solved. 

The  latter  type  of  resonance  was  discussed  in  detail  by  Poggio,  et 

al.  [12].  A  similar  phenomenon  in  acoustics  was  analyzed  by  Copley  [13] 

and  Schenck  [14].  Copley  showed  that  the  integral  representation  of  the 

velocity  potential  contained  the  product  of  1 / ( k  — k  )  and  a  surface 

m 

Integral,  where  k  is  the  free  space  wavenumber  and  k  an  eigen-value  of 

m 

the  interior  problem.  As  (k-k^)  approaches  zero,  the  surface  integral 
must  also  vanish  so  that  their  product  will  yield  a  finite  value  for  the 
velocity  potential.  Consequently,  it  becomes  increasingly  difficult  to 
calculate  accurately  the  ratio  of  two  vanishing  quantities  as  k 
approaches  k^.  Since  the  accuracy  of  the  computer  is  limited,  numerical 
inaccuracy  will  arise  at  the  resonant  frequency. 

In  electromagnetics,  a  similar  mechanism  was  noted  in  various 
integral  equation  formulations  [12,15,16,17].  Harrington  [16,17]  stated 
that  both  electric-field  and  magnetic-field  integral  equations  failed  at 
frequencies  corresponding  to  the  resonant  frequencies  of  the  interior 
problem.  The  present  algorithms  employ  the  magnetic  field  integral 
equation  and  are  therefore  expected  to  exhibit  difficulties  at  the 
interior  resonance  frequencies.  In  fact,  a  careful  search  near  the 
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first  eigen  frequency  of  a  sphere  of  radius  "a"  located  at  ka=2.744 
revealed  an  erratic  behavior  in  the  calculated  bac.kscatter  cross-section 
around  the  value  of  ka*2.89,  as  shown  in  Figure  4. 

The  distinctive  feature  of  the  internal  resonance  in  Figure  4  is 
the  rapidity  and  sharpness  of  the  variation  of  the  radar  cross-section 
versus  frequency.  Because  of  the  slight  frequency  shift  and  the 
narrowness  of  the  resonance  bandwidth,  this  resonance  escaped  an  earlier 
inquiry  into  this  phenomenon.  Although  the  resonance  has  shifted,  the 
radar  cross-section  curve  off  resonance  does  not  show  any  significant 
shift  one  way  or  the  other.  Therefore,  this  shift  of  resonance  can  not 
be  explained  as  being  due  to  the  effective  size  of  the  sphere 
simulation.  It  is  probably  due  to  the  surface-patch  algorithm,  which 
may  be  less  sensitive  to  resonance  phenomena  than  other  algorithms. 

It  is  now  clear  that  Internal  resonance  presents  a  numerical 
problem  in  the  present  algorithm  over  a  narrow  bandwidth  near  the 
resonance  frequency.  Correction  of  this  deficiency  can  be  carried  out 
by  following  either  of  the  two  approaches  outlined  by  Mautz  and 
Harrington  [16,17].  These  methods  involve  the  use  of  a  modified 
integral  equation,  the  combined  field  or  the  combined  source,  for  which 
the  solution  is  unique. 

F.  The  Basis  and  Weighting  Functions 

The  basis  function  used  in  the  present  approach  was  defined  in 
Equations  (3)  through  (7).  It  is  natural  to  ask  whether  other  basis 
functions  can  be  used  to  any  advantage.  For  example,  it  is  sometimes 
possible  to  orient  all  the  basis  vectors  perpendicular  to  a  fixed 
coordinate  axis  instead  of  the  incident  wave  polarization.  This  was 
tried  for  the  case  of  spheres  with  a  surprising  result.  The  output  data 
are  nearly  identical  whether  the  basis  vectors  are  oriented  according  to 
the  Incidence  polarization  or  a  fixed  coordinate.  No  other  geometry  was 
tested  because  of  the  difficulty  in  defining  a  basis  vector 
perpendicular  to  both  the  surface  normal  vector  and  the  fixed 
coordinate.  We  have  not  been  able  to  obtain  a  satisfactory  explanation 
for  this  discovery  during  our  brief  examination  of  this  problem. 
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Basis  functions  other  than  the  subsectional  pulse  function  have 
been  considered  in  the  present  study.  These  alternate  basis  functions 
included  the  triangular  and  sinusoidal  pyramid  functions.  The 
possibility  of  extending  these  subsectional  bases  over  more  than  one 
patch  was  also  considered.  Since  the  current  is  continuous  for  a 
conducting  body  having  a  closed  surface,  the  triangular  and  sinusoidal 
pyramid  bases  must  be  extended  to  cover  more  than  one  patch  in  order  to 
achieve  continuity  for  the  current. 

Figure  5  shows  a  triangular  pyramid  basis  function  peaked  at  the 
center  of  the  bisectors  of  the  center  patch.  Assuming  that  all  the 
patches  are  organized  to  share  one  side  with  each  of  the  three  adjacent 
triangles,  we  can  establish  a  basis  function  which  is  centered  in  each 
patch  and  declines  to  zero  at  the  centers  of  the  adjacent  patches.  Note 
that  the  shaded  and  unshaded  triangular  pyramids  in  Figure  5  only 
represent  the  magnitude  of  the  current,  not  its  polarization. 

One  problem  with  the  shaded  basis  function  from  the  shaded  to  the 
unshaded  pyramid  in  Figure  5  is  that  the  current  is  zero  on  the  vertices 
of  the  triangular  patch.  It  is  therefore  desirable  to  modify  the  basis 
function  by  extending  it  from  2,  3,  4,  to  2',  3',  and  4',  respectively. 
The  decay  of  the  basis  function  from  1  to  2,  3,  and  4,  (or  2',  3',  and 
4')  can  be  changed  from  linear  to  sinusoidal  if  it  offers  any 
computational  advantages. 

The  weighting  function  used  in  the  present  approach  is  the  Dirac- 
delta  function  according  to  Equation  (9).  If  a  flat  pulse  function  is 
used,  the  matrix  element  is  then  given  by 
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Equation  (41)  now  involves  a  double  integration.  The  integration  with 
respect  to  s'  can  be  approximately  expressed  in  closed  form  according  to 
Equations  (38)  -  (40).  The  next  step,  integration  with  respect  to  s, 
appears  difficult  and  must  be  computed  numerically. 

The  Dirac-delta  function  can  be  used  as  a  weighting  function  for 
the  triangular  pyramid  basis  function.  In  fact,  six  out  of  the  seven 
major  wire  analysis  programs  use  Dirac-delta  functions  as  the  weighting 
function  and  this  method  is  referred  to  as  collocation  (18].  In  this 
manner,  the  complexity  of  the  triangular  pyramid  basis  function  can  be 
drastically  reduced. 
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SECTION  III 
NUMERICAL  TESTING 


Based  on  the  study  described  in  the  preceeding  section,  two 
computer  programs,  SCT51  and  SCT52,  were  developed  for  numerical 
testing.  These  two  computer  programs  have  been  described  in  a  separate 
document  [19]  and  thus  their  descriptions  will  not  be  repeated  herein. 
The  program  SCT51  applies  to  arbitrarily-shaped  scatterers  without 
matrix  integration  and  the  program  SCT52  applies  to  symmetrical 
scatterers  with  the  incident  wave  propagating  along  the  plane  of 
symmetry. 

Numerical  computations  have  been  performed  previously  for 
conducting  spheres  of  various  sizes  [4],  Thus,  it  is  of  interest  to 
consider  other  geometries  and  the  scatterers  studied  in  this  report 
include  prolate  spheroids  and  finite  length  circular  cylinders  and 
rectangular  cylinders.  The  results  of  the  computations  were  then 
compared  with  existing  data  in  the  literature.  While  excellent 
agreement  between  the  data  and  the  computations  has  been  observed,  there 
was  an  apparent  discrepancy  with  respect  to  polarization  in  some  of  the 
cases  which  were  considered.  Namely,  our  TE  (Transverse  Electric) 
calculation  may  agree  with  known  data  of  TM  (Transverse  Magnetic)  case, 
and  our  TM  case  may  agree  with  the  known  TE  data  in  the  literature.  In 
fact,  our  results  agree  with  some  sources  and  disagree  with  some  other 
sources  as  far  as  the  incidence  polarization  is  concerned.  We  have 
examined  the  set-up  of  the  incidence  wave  in  our  computer  program  and 
have  not  found  any  error.  At  this  point,  we  tentatively  assume  that  the 
discrepancy  is  due  to  confusion  in  data  presentation  with  respect  to 
denoting  the  polarization.  A  further  examination  of  the  literature  may 
resolve  this  apparent  notation  problem. 

The  incidence  polarization  problem  emerged  only  when  it  affected 
the  results.  In  many  instances,  the  geometries  are  symmetrical  and  the 
scattering  characteristics  are  identical  in  the  two  principal  planes. 
There  is  also  a  possibility  of  error  in  the  definitions  of  the  axes  to 
which  the  TE  and  TM  modes  are  referred.  The  axis  was  defined  in  all  the 
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data  presented  in  the  literature  but  this  definition  could  be  a 
potential  source  of  error. 

In  order  to  resolve  the  polarization  problem,  a  series  of 
measurements  is  planned  in  early  1980  when  the  compact  range  at  Georgia 
Tech  will  be  available  for  scattering  measurements. 

With  the  foregoing  discussion  in  mind,  we  proceed  to  present  the 
numerical  results,  which  when  we  ignore  the  discrepancies  in 
polarization  in  some  cases,  are  in  good  agreement  with  existing  data  in 
the  literature. 

The  computation  for  prolate  spheroids  employed  a  54-point,  96-patch 
spheroidal  structure  described  in  Reference  4,  which  is  easily  elongated 
in  one  dimension  to  generate  the  desired  prolate  spheroid.  Figures  6 
through  8  show  the  comparison  of  the  back-scattering  cross-section 
resulting  from  the  present  calculation  and  the  data  from  Moffatt,  et  al. 
(20,21].  The  minor  axis  a  of  the  patched  prolate  spheroid  is  determined 
by  its  "effective"  value,  which  is  calculated  by  dividing  the 
circumference  by  2ir. 

The  backscatter  from  a  finite  length  circular  cylinder  was 
calculated  using  the  program  SCT52  and  the  results  are  compared  with 
data  in  the  literature  (22]  for  E-plane  and  H-plane  incidence.  Again 
the  question  of  effective  diameter  of  the  cylinder  arises.  Figure  9 
shows  how  the  finite  circular  cylinder  was  simulated  with  triangular 
patches.  Figure  10  shows  how  the  patched  cylinder  is  related  to  the 
circular  cylinder.  Obviously,  a  cross-section  Inscribed  on  circle  is  a 
less  reasonable  choice  than  the  cross-section  area  corrected  to  a 
circle.  Indeed,  better  results  were  obtained  by  the  method  of  cross- 
section  area  corrected  to  a  circle  as  shown  in  Fgures  11  and  12. 

Figures  13  and  14  show  the  comparison  of  the  computed  results  for  a 
finite  rectangular  cylinder  with  the  calculation  by  Tsai,  et  al  (23,24], 
who  referred  to  the  rectangular  cylinders  as  boxes.  SCT52,  which  takes 
advantage  of  the  symmetry  of  the  problem,  was  employed  in  this 
calculation.  Tsai's  calculation  also  utilized  the  symmetry  of  the 
scatterer.  Note  that  a  60  x  60  matrix  was  used  in  the  present 
calculation  while  a  100  x  100  matrix  was  used  by  Tsai. 


25 


*C3/i* 


Figure  7.  Rackscattering  cross  section  of  a  conducting 
prolate  spheroid  with  incidence  perpendicular 
to  the  major  axis  —  TE  polarization. 


CROSS-SECTION  AREA 
CORRECTED  TO  CIRCLE 


Figure  10.  Cross  sectional  views  showing  two  methods  in 

handling  the  effective  diameter  of  the  circular 
cylinder. 
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RADAR  CROSS-SECTION  (X*) 


Figure  14.  Comparison  of  the  computed  radar  cross  section  of  a 

conducting  rectangular  cylinder  with  other  known  data. 
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Several  salient  features  were  observed  In  the  numerical  testing  of 
the  Georgia  Tech  computer  program.  The  computer  program  is  capable  of 
generating  accurate  scattering  data  as  long  as  the  patch  sizes  are  such 
that  the  centers  of  the  adjacent  triangular  patches  are  spaced  a  quarter 
wavelength  or  less.  Geometries  involving  sharp  angles  such  as  the  boxes 
and  finite  circular  cylinders  can  be  handled  without  degradation  in 
accuracy.  The  computational  speed  and  central  memory  requirement  are 
dependent  on  the  number  of  patches  used  in  the  simulation.  For  a  96- 
patch  spheroid,  which  nearly  occupies  the  full  central  memory  of  the  CDC 
Cyber  74  computer  at  Georgia  Tech,  one  run  using  SCT51  at  a  single 
frequency  and  one  incidence  angle  takes  about  150  CPU  seconds.  When 
SCT52  is  used  for  a  symmetrical  geometry,  both  the  central  memory  and 
CPU  time  is  reduced  by  three-fourths. 

It  is  also  of  interest  to  compare  the  numerical  advantages  of  the 
surface-patch  modeling  and  the  wire-grid  modeling  approaches.  Although 
a  number  of  wire-grid  algorithms  have  been  developed  [18],  there  appears 
to  be  no  calculations  with  wire  algorithms  for  the  geometries  presented  in 
this  report.  A  wire-grid  model  for  a  sphere,  originally  calculated  by 
J.H.  Richmond,  was  reported  in  Reference  12.  This  wire-grid  model 
employed  an  earlier  algorithm  by  Richmond,  yielding  only  fair  results. 

A  more  critical  comparis  m  was  performed  at  Georgia  Tech  [4]  with 
Richmond's  latest  refined  wire  algorithm  [6].  It  was  shown  In  this 
comparison  that  the  present  surface-patch  approach  was  able  to  produce 
data  of  better  accuracy  while  using  fewer  linear  equations  and  less 
execution  time.  The  wire  approach  employed  194  equations,  and  required 
a  27-second  compilation  time  and  a  237-second  execution  time.  The 
surface-patch  approach  employed  192  equations,  and  required  a  8-second 
compilation  time  and  a  150-second  execution  time.  However,  Richmond's 
reaction  integral  equation  approach  has  one  major  advantage.  Its  matrix 
is  symmetrical  and  as  a  result  its  central  memory  requirement  is  only 
about  60  percent  of  that  for  the  surface-patch  algorithm.  This 
advantage  of  the  reaction  integral  equation  algorithm  is  apparently  not 
shared  by  other  wire  algorithms  and  was  not  mentioned  in  the  detailed 
comparison  study  of  various  wire-grid  approaches  [18], 
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SCT52,  the  surface-patch  algorithm  for  symmetrical  scatterers, 
reduces  the  execution  time  and  central  memory  requirements  by  75 
percent.  When  there  is  two-plane  symmetry,  this  reduction  is  94  percent. 
Thus,  for  the  sphere,  in  which  two-plane  symmetry  exists,  the  exectuion 
time  is  9  seconds  and  the  central  memory  size  will  be  one  sixteenth  of 
that  required  by  SCT51  or  one-eighth  of  that  required  by  Richmond's  wire 
algorithm. 

No  comparison  has  been  made  between  wire  algorithms  and  the  surface- 
patch  algorithms  for  geometries  other  than  the  sphere.  The  accuracy  and 
convergence  of  the  surface-patch  algorithm  discussed  herein  surpass 
those  of  the  wire  algorithms  so  that  a  detailed  comparison  with  wire 
algorithms  is  not  of  high  priority  in  the  research  program.  The  effort 
in  generating  correct  data  for  various  geometries  using  the  wire 
approach  is  estimated  at  one  to  two  person-months.  There  will  also  be 
the  question  of  what  size  of  diameter  to  choose  for  the  wires  and  how 
these  wires  should  be  organized  in  the  simulation.  The  fact  that  the 
scattering  characteristics  depend  on  the  wire  radius  to  be  chosen  has 
already  seriously  hindered  the  wire-grid  approach  for  surface  modeling. 
Now  that  the  surface-patch  approach  has  demonstrated  its  accuracy  for 
several  important  types  of  scatterers,  including  some  with  sharp  edges, 
the  advantages  of  the  surface-patch  approach  over  the  wire-gird 
approaches  are  evident. 
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SECTION  IV 

DEVELOPMENT  OF  A  COMPUTER  PROGRAM  FOR  THIN-SHELL  PROBLEMS 


It  is  well  known  that  the  magnetic  field  integral  equation 
encounters  numerical  instability  problems  when  used  for  thin-shell 
scatterers  [25],  This  difficulty  can  be  overcome  by  using  the  reaction 
integral  equation  [26]  and  electric  field  integral  equation  [27]  or 
others  [25,28,29],  However,  computations  in  the  literature  invariably 
are  limited  to  the  simple  case  of  a  rectangular  conducting  plate  and  a 
dihedral  corner  reflector.  There  is  no  indication  that  any  of  these 
methods  can  be  readily  applied  to  other  geometries. 

In  the  process  of  developing  a  computer  program  for  the  scattering 
of  arbitrarily-shaped  thin  shell  structures,  two  approaches  were  taken 
in  this  study.  A  reaction  integral  equation  approach  was  initially 
investigated  but  was  not  completed  because  of  large  CPU  time 
requirements  in  the  numerical  computation  of  matrix  elements.  An 
electric  field  integral  equation  approach  was  then  adopted  and  a 
computer  program  coded  and  debugged.  Both  of  these  approaches  are 
discussed  in  the  following  paragraphs. 

A,  Reaction  Integral  Equation  Approach 

The  reaction  integral  equation  has  been  applied  to  the  cases  of  a 
rectangular  plate  and  a  dihedral  corner  reflector  126],  The  possibility 
of  using  this  approach  for  arbitrarily-shaped  scatterers  is  examined 
herein.  The  major  difficulty  in  this  approach  lies  in  the  difficulties 
in  the  integration  to  obtain  the  matrix  elements  given  by 
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Strictly  speaking,  Equation  (42)  involves  triple  integrations  as 
E^  must  also  be  evaluated  by  an  integration  process.  Presently,  there 
appears  to  be  no  closed-form  expression  for  the  field  due  to  a 
triangular  current  patch.  Even  in  the  case  of  a  finite  line  source,  the 
sinusoidally  excited  dipole  is  probably  the  only  one  with  a  simple 
closed-fora  expression  for  the  near-zone  field  [6],  In  order  to  compute 
in  Equation  (42),  two  integration  methods  were  tried.  The  first 
method  employs  modern  numerical  integration  techniques,  and  the  second 
method  is  based  on  the  approximation  of  the  surface  current  with  several 
line  current  elements. 

In  the  numerical  integration  technique,  the  definite  integral  is 
expanded  into  a  finite  series  which  can  be  computed  numerically  [30], 
Specifically,  the  integration  of  a  function  over  a  triangular  area  can 
be  carried  out  with  a  64-point  formula  [30],  64  points  in  the  triangle 

are  pre-selected  according  to  a  simple  arithmatic  formula  and  the  values 
of  the  integrand  at  these  points  are  then  computed.  The  value  of  the 
integral  is  then  obtained  by  summing  up  the  product  of  these  64  sampled 
integrands  and  a  predetermined  weighting  function  of  simple  arithmatic 
fora. 

The  accuracy  of  the  numerical  integration  depends  on  how  rapidly 
and  how  frequently  the  value  of  the  integrand  varies  in  the  area  of 
integration.  The  64-point  algorithm  was  checked  with  several  known 
functions  and  the  accuracy  of  this  algorithm  was  quite  impressive.  For 
example,  Table  I  shows  the  comparison  between  the  exact  values  and  the 
results  of  numerical  integration  for  the  integral 

sinujxdydx.  (43) 

As  can  be  seen,  even  when  o>“10,  the  64-point  algorithm  is  highly 
accurate . 

Although  this  integration  algorithm  is  highly  accurate,  it  is 
inefficient.  The  total  subroutine  contains  only  17  lines  of  short  and 
simple  arithmetic  expressions,  as  well  as  an  8x8  "DO"  loop.  It  takes 
about  0.08  CPU  second  to  run  one  case  in  Table  I.  Consequently, 
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ACCURACY  OF  THE  64-POINT  NUMERICAL  INTEGRATION  ALGORITHM 


Integral  tested  = 


sin  wx  dy  dx. 


u 

Exact 

Numerical  integration 

0 

0 

0 

1 

0.1585290152 

0.1585290 

10 

0.1054402111 

0.1054402 

100 

0.0100506366 

0.05372281 

1000 

0.0009991731 

-0.008012371 
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considerable  restraint  must  be  exercised  in  applying  this  technique  to 
the  evaluation  of  matrix  elements  in  Equation  (42).  Since  we  were 
unable  to  obtain  a  closed  form  expression  for  in  Equation  (42) 
analytically,  we  chose  to  approximate  the  patch  current  with  current 
filaments  as  shown  in  Figure  15.  Obviously,  the  higher  the  number  of 
filaments  included  in  the  process,  the  more  accurate  will  be  the 
approximation . 

The  expressions  for  the  electric  field  radiated  from  a  line  current 

having  a  sinusoidal  distribution  were  well  documented  by  Schelkunoff  and 

Friis  [31],  Because  of  the  symmetry  of  a  straight  line  current,  the 

radiated  field  is  constant  around  the  axis  of  the  line  current.  The 

radiated  field  therefore  consists  of  two  components,  one  parallel  to  the 

current  and  one  perpendicular  to  the  current.  There  is  no  <t>  component  if 

a  cylindrical  coordinate  is  assigned  with  the  current  along  the  z-axis. 

A  check  of  the  formulas  for  E  and  E  showed  that  E„  was  an  exact 

P  z  P 

expression.  However,  the  expression  for  E^  was  not  exact  as  implied. 

In  Equations  (12)  through  (40)  of  Reference  31,  the  term 


1  II 
jue  3z' 


(44) 


was  omitted  in  the  derivation.  The  exact  expression  should  read 


dE 

z 


lit 

3z  1 


e2«t] 


I  (z ' ) dz 1 


(45) 


It  was  noted  that  Equation  (45)  without  the  term  of  Equation  (44) 
is  Pocklington' s  integral  equation  [32],  On  the  other  hand.  Equation 
(45)  was  employed  by  Kyle  (33],  Richmond  [34],  Harrington  [35],  and 
Thiele  [36]  in  their  thin-wire  computations.  In  the  present  study, 
Pocklington' s  integral  was  chosen  so  that  all  the  formulations  for  a 
line  current  source  were  based  on  Reference  31. 

Figure  16  shows  a  comparison  between  the  present  calculation  and 
the  calculated  data  in  Reference  26  for  the  mutual  impedance  between  two 
rectangular  dipoles.  The  present  calculation  employed  three  filament 
approximations  for  triangular  surface  patches.  Each  rectangular  patch 
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-  9.484  -  j  6.574 

U- 


Figure  16.  Comparison  of  computed  mutual  impedances  between 
two  rectangular  dipoles. 
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was  divided  diagonally  Into  two  triangular  patches.  The  agreement  with 
data  of  Wang,  Richmond,  and  Gilbreath  [26]  was  good  as  long  as  the 
dipoles  are  spaced  one  quarter  wavelength  away.  For  closely  spaced 
dipoles,  more  current  filaments  are  needed  in  the  approximation. 

This  combined  analytical-numerical  method  made  it  possible  to 
reduce  drastically  the  computational  time  required  for  the  matrix 
element  to  about  1/20  of  the  time  needed  for  the  previous  numerical 
integration.  However,  the  requirement  for  computational  time  was  still 
prohibitively  high.  There  appears  to  be  no  readily  available  technique 
to  reduce  the  CPU  time  to  an  acceptable  level.  The  course  of  research 
was  then  turned  to  the  electric  field  integral  equation  approach. 

B,  Electric  Field  Integral  Equation  Approach 

The  electric  field  integral  equation  usually  takes  the  following 
form  [12] 

n  x  Emc(r)  =  n  x  £{-u>  (r  * ) <f>(r , r ' ) 

+  [V'  .  J  (r ' )  ] V' 4>(r  ,r ' ) }  ds'  (46) 

S  “S  - 

where  roost  of  notations  have  been  defined  in  Equation  (1)  and 

in  c 

E  (r)  «  incident  electric  field, 

v'  -l  +  t  _i_ 

s  L1  3t:  c2  3t2  ’ 

t^  ,  t^  are  two  orthonormal  unit  vectors  on  the  surface  S, 


We  denote  the  right-hand  side  of  Equation  (46)  with  an  operator  form 
X(J  (r'))  and  rewrite  Equation  (46)  as 


x(4d')) 


n  x 


(47) 
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We  then  expand  the  surface  current  J  (r1)  according  to  Equations 
(3)  through  (7)  and  employ  the  method  of  moments  [35] by  performing  the 
scalar  product  according  to  Equations  (9)  and  (10),  The  resulting 
system  of  linear  equations  is 
2  N 

£  E 

j=l  n= 1 


Ij  Aij 
n  mn 


‘  .  inc,  ,  „i 

=  n  x  E  <  r  )  *  U 
m  —  —  m  m 


(48) 


i=l ,  2 ;  m=l , .  . . N , 


where  for  m^n 


Alj  =  Aij 


lmn 


+  A 


ij 

2mn 


(49) 
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for  m=n 


(52) 
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In  the  foregoing  derivation,  the  surface  gradient  operator  7'  is 
approximated  by 


V'  J  (r  '  ) 
s  — n  — 


3  1 1 


at1 

n 


1  _l1  i 

n  n  1 

r  .  U1 

— n  n  1  — n 


(53) 


where  nl  and  n2  are  the  indices  for  two  patches  adjacent  to  n  and  are 

1  2 

chosen  for  the  evaluation  of  the  gradient  of  1  and  I  .  Note  that  the 

n  n 

accuracy  of  Equation  (53)  depends  on  the  local  curvature  of  the  surface 

S.  For  a  flat  plane,  Equation  (13)  is  most  accurate  because  the 

polarizations  of  the  first  basis  vectors  U  *  and  U  , *  are  identical;  so 

— n  — nl  2  2 

are  the  polarizations  for  the  second  basis  vectors  U  and  U  ,  .  These 
r  — n  — nl 

relations  also  hold  for  patch  n2.  r  .  and  r  „  in  Equation  (53)  are 

— nml  — nm2 

the  vectorial  distances  between  the  centers  of  patch  n  and  patches  nl 
and  n2,  respectively,  as  shown  in  Figure  17, 

When  one  side  of  the  n1"^  triangular  patch  is  an  edge,  the  basis 

vectors  on  this  patch  are  organized  to  orient  U  *  to  be  perpendicular 

1 

to  the  edge,  I  ,  being  perpendicular  to  the  edge,  is  therefore  zero  in 
— n 

magnitude.  r  ,  is  reduced  to  the  case  in  which  the  area  of  patch  nl  is 
— nnl 

zero.  In  other  words,  JLnni  runs  from  the  bisecting  piont  on  the  edge 
side  to  the  center  of  the  triangle. 

A  generalized  computer  program  based  on  the  foregoing  analysis  has 
been  coded.  There  are  many  similarities  between  the  new  thin-shell 
program  and  SCT51  program  using  the  magnetic  field  integral  equation. 
Numerical  testing  of  the  thin-shell  computer  program  will  be  conducted 
and  the  results  will  be  included  in  the  final  report. 
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SECTION  V 
CONCLUSIONS 


Numerical  testing  for  prolate  spheroids,  and  finite  circular  and 
rectangular  cylinders  has  demonstrated  that  the  surface  patch  algorithms 
are  efficient  and  accurate.  These  algorithms  are  also  capable  of 
handling  sharp  edges  in  rectangular  box  and  finite  cylinders.  As 
compared  with  the  wire-grid  algorithms,  the  surface  patch  algorithm  is 
more  efficient  with  respect  to  computation  speed  and  exhibits  improved 
numerical  accuracy,  but  has  the  disadvantage  of  unsymmetrical  matrices. 
Improvements  and  limitations  of  this  algorithm  have  also  been  discussed. 
Research  efforts  on  the  development  of  a  thin-shell  surface-patch 
algorithm  is  also  presented.  This  thin-shell  algorithm  is  being  tested 
numerically  for  its  accuracy  and  convergence. 
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MISSION  I 

°f 

Ram  Air  Development  Center  , 

RADC  plans  and  execute*  research,  development,  tut  and  a 

*  elected  acquisition  programs  l n  support  of  Command,  Control  C 
Communication*  and  Intelligence  (C3I)  activities.  Technical  , 


communication* ,  electromagnetic  guidance  and  control,  sur-  a 
veillance  oi  ground  and  aerospace  objects,  intelligence  data  % 
collection  and  handling,  information  system  technology,  C 

ionospheric  propagation,  solid  state  sciences,  microwave  3 

physics  and  electronic  reliability,  maintainability  and  0 

compatibility.  8 
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